// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// Copyright (C) 2012 Alexey Korepanov <kaikaikai@yandex.ru>
//
// This Source Code Form is subject to the terms of the Mozilla
// Public License v. 2.0. If a copy of the MPL was not distributed
// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.

#ifndef EIGEN_REAL_QZ_H
#define EIGEN_REAL_QZ_H

namespace Eigen {

/** \eigenvalues_module \ingroup Eigenvalues_Module
 *
 *
 * \class RealQZ
 *
 * \brief Performs a real QZ decomposition of a pair of square matrices
 *
 * \tparam _MatrixType the type of the matrix of which we are computing the
 * real QZ decomposition; this is expected to be an instantiation of the
 * Matrix class template.
 *
 * Given a real square matrices A and B, this class computes the real QZ
 * decomposition: \f$ A = Q S Z \f$, \f$ B = Q T Z \f$ where Q and Z are
 * real orthogonal matrixes, T is upper-triangular matrix, and S is upper
 * quasi-triangular matrix. An orthogonal matrix is a matrix whose
 * inverse is equal to its transpose, \f$ U^{-1} = U^T \f$. A quasi-triangular
 * matrix is a block-triangular matrix whose diagonal consists of 1-by-1
 * blocks and 2-by-2 blocks where further reduction is impossible due to
 * complex eigenvalues.
 *
 * The eigenvalues of the pencil \f$ A - z B \f$ can be obtained from
 * 1x1 and 2x2 blocks on the diagonals of S and T.
 *
 * Call the function compute() to compute the real QZ decomposition of a
 * given pair of matrices. Alternatively, you can use the
 * RealQZ(const MatrixType& B, const MatrixType& B, bool computeQZ)
 * constructor which computes the real QZ decomposition at construction
 * time. Once the decomposition is computed, you can use the matrixS(),
 * matrixT(), matrixQ() and matrixZ() functions to retrieve the matrices
 * S, T, Q and Z in the decomposition. If computeQZ==false, some time
 * is saved by not computing matrices Q and Z.
 *
 * Example: \include RealQZ_compute.cpp
 * Output: \include RealQZ_compute.out
 *
 * \note The implementation is based on the algorithm in "Matrix Computations"
 * by Gene H. Golub and Charles F. Van Loan, and a paper "An algorithm for
 * generalized eigenvalue problems" by C.B.Moler and G.W.Stewart.
 *
 * \sa class RealSchur, class ComplexSchur, class EigenSolver, class ComplexEigenSolver
 */

template<typename _MatrixType>
class RealQZ
{
  public:
	typedef _MatrixType MatrixType;
	enum
	{
		RowsAtCompileTime = MatrixType::RowsAtCompileTime,
		ColsAtCompileTime = MatrixType::ColsAtCompileTime,
		Options = MatrixType::Options,
		MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
		MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
	};
	typedef typename MatrixType::Scalar Scalar;
	typedef std::complex<typename NumTraits<Scalar>::Real> ComplexScalar;
	typedef Eigen::Index Index; ///< \deprecated since Eigen 3.3

	typedef Matrix<ComplexScalar, ColsAtCompileTime, 1, Options & ~RowMajor, MaxColsAtCompileTime, 1> EigenvalueType;
	typedef Matrix<Scalar, ColsAtCompileTime, 1, Options & ~RowMajor, MaxColsAtCompileTime, 1> ColumnVectorType;

	/** \brief Default constructor.
	 *
	 * \param [in] size  Positive integer, size of the matrix whose QZ decomposition will be computed.
	 *
	 * The default constructor is useful in cases in which the user intends to
	 * perform decompositions via compute().  The \p size parameter is only
	 * used as a hint. It is not an error to give a wrong \p size, but it may
	 * impair performance.
	 *
	 * \sa compute() for an example.
	 */
	explicit RealQZ(Index size = RowsAtCompileTime == Dynamic ? 1 : RowsAtCompileTime)
		: m_S(size, size)
		, m_T(size, size)
		, m_Q(size, size)
		, m_Z(size, size)
		, m_workspace(size * 2)
		, m_maxIters(400)
		, m_isInitialized(false)
		, m_computeQZ(true)
	{
	}

	/** \brief Constructor; computes real QZ decomposition of given matrices
	 *
	 * \param[in]  A          Matrix A.
	 * \param[in]  B          Matrix B.
	 * \param[in]  computeQZ  If false, A and Z are not computed.
	 *
	 * This constructor calls compute() to compute the QZ decomposition.
	 */
	RealQZ(const MatrixType& A, const MatrixType& B, bool computeQZ = true)
		: m_S(A.rows(), A.cols())
		, m_T(A.rows(), A.cols())
		, m_Q(A.rows(), A.cols())
		, m_Z(A.rows(), A.cols())
		, m_workspace(A.rows() * 2)
		, m_maxIters(400)
		, m_isInitialized(false)
		, m_computeQZ(true)
	{
		compute(A, B, computeQZ);
	}

	/** \brief Returns matrix Q in the QZ decomposition.
	 *
	 * \returns A const reference to the matrix Q.
	 */
	const MatrixType& matrixQ() const
	{
		eigen_assert(m_isInitialized && "RealQZ is not initialized.");
		eigen_assert(m_computeQZ && "The matrices Q and Z have not been computed during the QZ decomposition.");
		return m_Q;
	}

	/** \brief Returns matrix Z in the QZ decomposition.
	 *
	 * \returns A const reference to the matrix Z.
	 */
	const MatrixType& matrixZ() const
	{
		eigen_assert(m_isInitialized && "RealQZ is not initialized.");
		eigen_assert(m_computeQZ && "The matrices Q and Z have not been computed during the QZ decomposition.");
		return m_Z;
	}

	/** \brief Returns matrix S in the QZ decomposition.
	 *
	 * \returns A const reference to the matrix S.
	 */
	const MatrixType& matrixS() const
	{
		eigen_assert(m_isInitialized && "RealQZ is not initialized.");
		return m_S;
	}

	/** \brief Returns matrix S in the QZ decomposition.
	 *
	 * \returns A const reference to the matrix S.
	 */
	const MatrixType& matrixT() const
	{
		eigen_assert(m_isInitialized && "RealQZ is not initialized.");
		return m_T;
	}

	/** \brief Computes QZ decomposition of given matrix.
	 *
	 * \param[in]  A          Matrix A.
	 * \param[in]  B          Matrix B.
	 * \param[in]  computeQZ  If false, A and Z are not computed.
	 * \returns    Reference to \c *this
	 */
	RealQZ& compute(const MatrixType& A, const MatrixType& B, bool computeQZ = true);

	/** \brief Reports whether previous computation was successful.
	 *
	 * \returns \c Success if computation was successful, \c NoConvergence otherwise.
	 */
	ComputationInfo info() const
	{
		eigen_assert(m_isInitialized && "RealQZ is not initialized.");
		return m_info;
	}

	/** \brief Returns number of performed QR-like iterations.
	 */
	Index iterations() const
	{
		eigen_assert(m_isInitialized && "RealQZ is not initialized.");
		return m_global_iter;
	}

	/** Sets the maximal number of iterations allowed to converge to one eigenvalue
	 * or decouple the problem.
	 */
	RealQZ& setMaxIterations(Index maxIters)
	{
		m_maxIters = maxIters;
		return *this;
	}

  private:
	MatrixType m_S, m_T, m_Q, m_Z;
	Matrix<Scalar, Dynamic, 1> m_workspace;
	ComputationInfo m_info;
	Index m_maxIters;
	bool m_isInitialized;
	bool m_computeQZ;
	Scalar m_normOfT, m_normOfS;
	Index m_global_iter;

	typedef Matrix<Scalar, 3, 1> Vector3s;
	typedef Matrix<Scalar, 2, 1> Vector2s;
	typedef Matrix<Scalar, 2, 2> Matrix2s;
	typedef JacobiRotation<Scalar> JRs;

	void hessenbergTriangular();
	void computeNorms();
	Index findSmallSubdiagEntry(Index iu);
	Index findSmallDiagEntry(Index f, Index l);
	void splitOffTwoRows(Index i);
	void pushDownZero(Index z, Index f, Index l);
	void step(Index f, Index l, Index iter);

}; // RealQZ

/** \internal Reduces S and T to upper Hessenberg - triangular form */
template<typename MatrixType>
void
RealQZ<MatrixType>::hessenbergTriangular()
{

	const Index dim = m_S.cols();

	// perform QR decomposition of T, overwrite T with R, save Q
	HouseholderQR<MatrixType> qrT(m_T);
	m_T = qrT.matrixQR();
	m_T.template triangularView<StrictlyLower>().setZero();
	m_Q = qrT.householderQ();
	// overwrite S with Q* S
	m_S.applyOnTheLeft(m_Q.adjoint());
	// init Z as Identity
	if (m_computeQZ)
		m_Z = MatrixType::Identity(dim, dim);
	// reduce S to upper Hessenberg with Givens rotations
	for (Index j = 0; j <= dim - 3; j++) {
		for (Index i = dim - 1; i >= j + 2; i--) {
			JRs G;
			// kill S(i,j)
			if (m_S.coeff(i, j) != 0) {
				G.makeGivens(m_S.coeff(i - 1, j), m_S.coeff(i, j), &m_S.coeffRef(i - 1, j));
				m_S.coeffRef(i, j) = Scalar(0.0);
				m_S.rightCols(dim - j - 1).applyOnTheLeft(i - 1, i, G.adjoint());
				m_T.rightCols(dim - i + 1).applyOnTheLeft(i - 1, i, G.adjoint());
				// update Q
				if (m_computeQZ)
					m_Q.applyOnTheRight(i - 1, i, G);
			}
			// kill T(i,i-1)
			if (m_T.coeff(i, i - 1) != Scalar(0)) {
				G.makeGivens(m_T.coeff(i, i), m_T.coeff(i, i - 1), &m_T.coeffRef(i, i));
				m_T.coeffRef(i, i - 1) = Scalar(0.0);
				m_S.applyOnTheRight(i, i - 1, G);
				m_T.topRows(i).applyOnTheRight(i, i - 1, G);
				// update Z
				if (m_computeQZ)
					m_Z.applyOnTheLeft(i, i - 1, G.adjoint());
			}
		}
	}
}

/** \internal Computes vector L1 norms of S and T when in Hessenberg-Triangular form already */
template<typename MatrixType>
inline void
RealQZ<MatrixType>::computeNorms()
{
	const Index size = m_S.cols();
	m_normOfS = Scalar(0.0);
	m_normOfT = Scalar(0.0);
	for (Index j = 0; j < size; ++j) {
		m_normOfS += m_S.col(j).segment(0, (std::min)(size, j + 2)).cwiseAbs().sum();
		m_normOfT += m_T.row(j).segment(j, size - j).cwiseAbs().sum();
	}
}

/** \internal Look for single small sub-diagonal element S(res, res-1) and return res (or 0) */
template<typename MatrixType>
inline Index
RealQZ<MatrixType>::findSmallSubdiagEntry(Index iu)
{
	using std::abs;
	Index res = iu;
	while (res > 0) {
		Scalar s = abs(m_S.coeff(res - 1, res - 1)) + abs(m_S.coeff(res, res));
		if (s == Scalar(0.0))
			s = m_normOfS;
		if (abs(m_S.coeff(res, res - 1)) < NumTraits<Scalar>::epsilon() * s)
			break;
		res--;
	}
	return res;
}

/** \internal Look for single small diagonal element T(res, res) for res between f and l, and return res (or f-1)  */
template<typename MatrixType>
inline Index
RealQZ<MatrixType>::findSmallDiagEntry(Index f, Index l)
{
	using std::abs;
	Index res = l;
	while (res >= f) {
		if (abs(m_T.coeff(res, res)) <= NumTraits<Scalar>::epsilon() * m_normOfT)
			break;
		res--;
	}
	return res;
}

/** \internal decouple 2x2 diagonal block in rows i, i+1 if eigenvalues are real */
template<typename MatrixType>
inline void
RealQZ<MatrixType>::splitOffTwoRows(Index i)
{
	using std::abs;
	using std::sqrt;
	const Index dim = m_S.cols();
	if (abs(m_S.coeff(i + 1, i)) == Scalar(0))
		return;
	Index j = findSmallDiagEntry(i, i + 1);
	if (j == i - 1) {
		// block of (S T^{-1})
		Matrix2s STi = m_T.template block<2, 2>(i, i).template triangularView<Upper>().template solve<OnTheRight>(
			m_S.template block<2, 2>(i, i));
		Scalar p = Scalar(0.5) * (STi(0, 0) - STi(1, 1));
		Scalar q = p * p + STi(1, 0) * STi(0, 1);
		if (q >= 0) {
			Scalar z = sqrt(q);
			// one QR-like iteration for ABi - lambda I
			// is enough - when we know exact eigenvalue in advance,
			// convergence is immediate
			JRs G;
			if (p >= 0)
				G.makeGivens(p + z, STi(1, 0));
			else
				G.makeGivens(p - z, STi(1, 0));
			m_S.rightCols(dim - i).applyOnTheLeft(i, i + 1, G.adjoint());
			m_T.rightCols(dim - i).applyOnTheLeft(i, i + 1, G.adjoint());
			// update Q
			if (m_computeQZ)
				m_Q.applyOnTheRight(i, i + 1, G);

			G.makeGivens(m_T.coeff(i + 1, i + 1), m_T.coeff(i + 1, i));
			m_S.topRows(i + 2).applyOnTheRight(i + 1, i, G);
			m_T.topRows(i + 2).applyOnTheRight(i + 1, i, G);
			// update Z
			if (m_computeQZ)
				m_Z.applyOnTheLeft(i + 1, i, G.adjoint());

			m_S.coeffRef(i + 1, i) = Scalar(0.0);
			m_T.coeffRef(i + 1, i) = Scalar(0.0);
		}
	} else {
		pushDownZero(j, i, i + 1);
	}
}

/** \internal use zero in T(z,z) to zero S(l,l-1), working in block f..l */
template<typename MatrixType>
inline void
RealQZ<MatrixType>::pushDownZero(Index z, Index f, Index l)
{
	JRs G;
	const Index dim = m_S.cols();
	for (Index zz = z; zz < l; zz++) {
		// push 0 down
		Index firstColS = zz > f ? (zz - 1) : zz;
		G.makeGivens(m_T.coeff(zz, zz + 1), m_T.coeff(zz + 1, zz + 1));
		m_S.rightCols(dim - firstColS).applyOnTheLeft(zz, zz + 1, G.adjoint());
		m_T.rightCols(dim - zz).applyOnTheLeft(zz, zz + 1, G.adjoint());
		m_T.coeffRef(zz + 1, zz + 1) = Scalar(0.0);
		// update Q
		if (m_computeQZ)
			m_Q.applyOnTheRight(zz, zz + 1, G);
		// kill S(zz+1, zz-1)
		if (zz > f) {
			G.makeGivens(m_S.coeff(zz + 1, zz), m_S.coeff(zz + 1, zz - 1));
			m_S.topRows(zz + 2).applyOnTheRight(zz, zz - 1, G);
			m_T.topRows(zz + 1).applyOnTheRight(zz, zz - 1, G);
			m_S.coeffRef(zz + 1, zz - 1) = Scalar(0.0);
			// update Z
			if (m_computeQZ)
				m_Z.applyOnTheLeft(zz, zz - 1, G.adjoint());
		}
	}
	// finally kill S(l,l-1)
	G.makeGivens(m_S.coeff(l, l), m_S.coeff(l, l - 1));
	m_S.applyOnTheRight(l, l - 1, G);
	m_T.applyOnTheRight(l, l - 1, G);
	m_S.coeffRef(l, l - 1) = Scalar(0.0);
	// update Z
	if (m_computeQZ)
		m_Z.applyOnTheLeft(l, l - 1, G.adjoint());
}

/** \internal QR-like iterative step for block f..l */
template<typename MatrixType>
inline void
RealQZ<MatrixType>::step(Index f, Index l, Index iter)
{
	using std::abs;
	const Index dim = m_S.cols();

	// x, y, z
	Scalar x, y, z;
	if (iter == 10) {
		// Wilkinson ad hoc shift
		const Scalar a11 = m_S.coeff(f + 0, f + 0), a12 = m_S.coeff(f + 0, f + 1), a21 = m_S.coeff(f + 1, f + 0),
					 a22 = m_S.coeff(f + 1, f + 1), a32 = m_S.coeff(f + 2, f + 1), b12 = m_T.coeff(f + 0, f + 1),
					 b11i = Scalar(1.0) / m_T.coeff(f + 0, f + 0), b22i = Scalar(1.0) / m_T.coeff(f + 1, f + 1),
					 a87 = m_S.coeff(l - 1, l - 2), a98 = m_S.coeff(l - 0, l - 1),
					 b77i = Scalar(1.0) / m_T.coeff(l - 2, l - 2), b88i = Scalar(1.0) / m_T.coeff(l - 1, l - 1);
		Scalar ss = abs(a87 * b77i) + abs(a98 * b88i), lpl = Scalar(1.5) * ss, ll = ss * ss;
		x = ll + a11 * a11 * b11i * b11i - lpl * a11 * b11i + a12 * a21 * b11i * b22i -
			a11 * a21 * b12 * b11i * b11i * b22i;
		y = a11 * a21 * b11i * b11i - lpl * a21 * b11i + a21 * a22 * b11i * b22i - a21 * a21 * b12 * b11i * b11i * b22i;
		z = a21 * a32 * b11i * b22i;
	} else if (iter == 16) {
		// another exceptional shift
		x = m_S.coeff(f, f) / m_T.coeff(f, f) - m_S.coeff(l, l) / m_T.coeff(l, l) +
			m_S.coeff(l, l - 1) * m_T.coeff(l - 1, l) / (m_T.coeff(l - 1, l - 1) * m_T.coeff(l, l));
		y = m_S.coeff(f + 1, f) / m_T.coeff(f, f);
		z = 0;
	} else if (iter > 23 && !(iter % 8)) {
		// extremely exceptional shift
		x = internal::random<Scalar>(-1.0, 1.0);
		y = internal::random<Scalar>(-1.0, 1.0);
		z = internal::random<Scalar>(-1.0, 1.0);
	} else {
		// Compute the shifts: (x,y,z,0...) = (AB^-1 - l1 I) (AB^-1 - l2 I) e1
		// where l1 and l2 are the eigenvalues of the 2x2 matrix C = U V^-1 where
		// U and V are 2x2 bottom right sub matrices of A and B. Thus:
		//  = AB^-1AB^-1 + l1 l2 I - (l1+l2)(AB^-1)
		//  = AB^-1AB^-1 + det(M) - tr(M)(AB^-1)
		// Since we are only interested in having x, y, z with a correct ratio, we have:
		const Scalar a11 = m_S.coeff(f, f), a12 = m_S.coeff(f, f + 1), a21 = m_S.coeff(f + 1, f),
					 a22 = m_S.coeff(f + 1, f + 1), a32 = m_S.coeff(f + 2, f + 1),

					 a88 = m_S.coeff(l - 1, l - 1), a89 = m_S.coeff(l - 1, l), a98 = m_S.coeff(l, l - 1),
					 a99 = m_S.coeff(l, l),

					 b11 = m_T.coeff(f, f), b12 = m_T.coeff(f, f + 1), b22 = m_T.coeff(f + 1, f + 1),

					 b88 = m_T.coeff(l - 1, l - 1), b89 = m_T.coeff(l - 1, l), b99 = m_T.coeff(l, l);

		x = ((a88 / b88 - a11 / b11) * (a99 / b99 - a11 / b11) - (a89 / b99) * (a98 / b88) +
			 (a98 / b88) * (b89 / b99) * (a11 / b11)) *
				(b11 / a21) +
			a12 / b22 - (a11 / b11) * (b12 / b22);
		y = (a22 / b22 - a11 / b11) - (a21 / b11) * (b12 / b22) - (a88 / b88 - a11 / b11) - (a99 / b99 - a11 / b11) +
			(a98 / b88) * (b89 / b99);
		z = a32 / b22;
	}

	JRs G;

	for (Index k = f; k <= l - 2; k++) {
		// variables for Householder reflections
		Vector2s essential2;
		Scalar tau, beta;

		Vector3s hr(x, y, z);

		// Q_k to annihilate S(k+1,k-1) and S(k+2,k-1)
		hr.makeHouseholderInPlace(tau, beta);
		essential2 = hr.template bottomRows<2>();
		Index fc = (std::max)(k - 1, Index(0)); // first col to update
		m_S.template middleRows<3>(k).rightCols(dim - fc).applyHouseholderOnTheLeft(
			essential2, tau, m_workspace.data());
		m_T.template middleRows<3>(k).rightCols(dim - fc).applyHouseholderOnTheLeft(
			essential2, tau, m_workspace.data());
		if (m_computeQZ)
			m_Q.template middleCols<3>(k).applyHouseholderOnTheRight(essential2, tau, m_workspace.data());
		if (k > f)
			m_S.coeffRef(k + 2, k - 1) = m_S.coeffRef(k + 1, k - 1) = Scalar(0.0);

		// Z_{k1} to annihilate T(k+2,k+1) and T(k+2,k)
		hr << m_T.coeff(k + 2, k + 2), m_T.coeff(k + 2, k), m_T.coeff(k + 2, k + 1);
		hr.makeHouseholderInPlace(tau, beta);
		essential2 = hr.template bottomRows<2>();
		{
			Index lr = (std::min)(k + 4, dim); // last row to update
			Map<Matrix<Scalar, Dynamic, 1>> tmp(m_workspace.data(), lr);
			// S
			tmp = m_S.template middleCols<2>(k).topRows(lr) * essential2;
			tmp += m_S.col(k + 2).head(lr);
			m_S.col(k + 2).head(lr) -= tau * tmp;
			m_S.template middleCols<2>(k).topRows(lr) -= (tau * tmp) * essential2.adjoint();
			// T
			tmp = m_T.template middleCols<2>(k).topRows(lr) * essential2;
			tmp += m_T.col(k + 2).head(lr);
			m_T.col(k + 2).head(lr) -= tau * tmp;
			m_T.template middleCols<2>(k).topRows(lr) -= (tau * tmp) * essential2.adjoint();
		}
		if (m_computeQZ) {
			// Z
			Map<Matrix<Scalar, 1, Dynamic>> tmp(m_workspace.data(), dim);
			tmp = essential2.adjoint() * (m_Z.template middleRows<2>(k));
			tmp += m_Z.row(k + 2);
			m_Z.row(k + 2) -= tau * tmp;
			m_Z.template middleRows<2>(k) -= essential2 * (tau * tmp);
		}
		m_T.coeffRef(k + 2, k) = m_T.coeffRef(k + 2, k + 1) = Scalar(0.0);

		// Z_{k2} to annihilate T(k+1,k)
		G.makeGivens(m_T.coeff(k + 1, k + 1), m_T.coeff(k + 1, k));
		m_S.applyOnTheRight(k + 1, k, G);
		m_T.applyOnTheRight(k + 1, k, G);
		// update Z
		if (m_computeQZ)
			m_Z.applyOnTheLeft(k + 1, k, G.adjoint());
		m_T.coeffRef(k + 1, k) = Scalar(0.0);

		// update x,y,z
		x = m_S.coeff(k + 1, k);
		y = m_S.coeff(k + 2, k);
		if (k < l - 2)
			z = m_S.coeff(k + 3, k);
	} // loop over k

	// Q_{n-1} to annihilate y = S(l,l-2)
	G.makeGivens(x, y);
	m_S.applyOnTheLeft(l - 1, l, G.adjoint());
	m_T.applyOnTheLeft(l - 1, l, G.adjoint());
	if (m_computeQZ)
		m_Q.applyOnTheRight(l - 1, l, G);
	m_S.coeffRef(l, l - 2) = Scalar(0.0);

	// Z_{n-1} to annihilate T(l,l-1)
	G.makeGivens(m_T.coeff(l, l), m_T.coeff(l, l - 1));
	m_S.applyOnTheRight(l, l - 1, G);
	m_T.applyOnTheRight(l, l - 1, G);
	if (m_computeQZ)
		m_Z.applyOnTheLeft(l, l - 1, G.adjoint());
	m_T.coeffRef(l, l - 1) = Scalar(0.0);
}

template<typename MatrixType>
RealQZ<MatrixType>&
RealQZ<MatrixType>::compute(const MatrixType& A_in, const MatrixType& B_in, bool computeQZ)
{

	const Index dim = A_in.cols();

	eigen_assert(A_in.rows() == dim && A_in.cols() == dim && B_in.rows() == dim && B_in.cols() == dim &&
				 "Need square matrices of the same dimension");

	m_isInitialized = true;
	m_computeQZ = computeQZ;
	m_S = A_in;
	m_T = B_in;
	m_workspace.resize(dim * 2);
	m_global_iter = 0;

	// entrance point: hessenberg triangular decomposition
	hessenbergTriangular();
	// compute L1 vector norms of T, S into m_normOfS, m_normOfT
	computeNorms();

	Index l = dim - 1, f, local_iter = 0;

	while (l > 0 && local_iter < m_maxIters) {
		f = findSmallSubdiagEntry(l);
		// now rows and columns f..l (including) decouple from the rest of the problem
		if (f > 0)
			m_S.coeffRef(f, f - 1) = Scalar(0.0);
		if (f == l) // One root found
		{
			l--;
			local_iter = 0;
		} else if (f == l - 1) // Two roots found
		{
			splitOffTwoRows(f);
			l -= 2;
			local_iter = 0;
		} else // No convergence yet
		{
			// if there's zero on diagonal of T, we can isolate an eigenvalue with Givens rotations
			Index z = findSmallDiagEntry(f, l);
			if (z >= f) {
				// zero found
				pushDownZero(z, f, l);
			} else {
				// We are sure now that S.block(f,f, l-f+1,l-f+1) is underuced upper-Hessenberg
				// and T.block(f,f, l-f+1,l-f+1) is invertible uper-triangular, which allows to
				// apply a QR-like iteration to rows and columns f..l.
				step(f, l, local_iter);
				local_iter++;
				m_global_iter++;
			}
		}
	}
	// check if we converged before reaching iterations limit
	m_info = (local_iter < m_maxIters) ? Success : NoConvergence;

	// For each non triangular 2x2 diagonal block of S,
	//    reduce the respective 2x2 diagonal block of T to positive diagonal form using 2x2 SVD.
	// This step is not mandatory for QZ, but it does help further extraction of eigenvalues/eigenvectors,
	// and is in par with Lapack/Matlab QZ.
	if (m_info == Success) {
		for (Index i = 0; i < dim - 1; ++i) {
			if (m_S.coeff(i + 1, i) != Scalar(0)) {
				JacobiRotation<Scalar> j_left, j_right;
				internal::real_2x2_jacobi_svd(m_T, i, i + 1, &j_left, &j_right);

				// Apply resulting Jacobi rotations
				m_S.applyOnTheLeft(i, i + 1, j_left);
				m_S.applyOnTheRight(i, i + 1, j_right);
				m_T.applyOnTheLeft(i, i + 1, j_left);
				m_T.applyOnTheRight(i, i + 1, j_right);
				m_T(i + 1, i) = m_T(i, i + 1) = Scalar(0);

				if (m_computeQZ) {
					m_Q.applyOnTheRight(i, i + 1, j_left.transpose());
					m_Z.applyOnTheLeft(i, i + 1, j_right.transpose());
				}

				i++;
			}
		}
	}

	return *this;
} // end compute

} // end namespace Eigen

#endif // EIGEN_REAL_QZ
